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Abstract. An intriguing and unexpected result for students learning numerical analysis is 
that Newton's method, applied to the simple polynomial z 3 — 1 = in the complex plane, leads to 
intricately interwoven basins of attraction of the roots. As an example of an interesting open question 
that may help to stimulate student interest in numerical analysis, we investigate the question of 
whether a damping method, which is designed to increase the likelihood of convergence for Newton's 
method, modifies the fractal structure of the basin boundaries. The overlap of the frontiers of 
numerical analysis and nonlinear dynamics provides many other problems that can help to make 
numerical analysis courses interesting. 
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1. Introduction. Numerical analysis is an immensely exciting area of research 
with many intellectual and practical benefits and yet this excitement sometimes 
doesn't come across in courses and textbooks. One reason is that courses and texts 
sometimes do not indicate some of the frontiers of the field and why these frontiers are 
interesting. Perhaps more importantly, few introductory courses and texts make stu- 
dents aware that they are actually close to, and can contribute to, research frontiers 
even while participating in an introductory course. The ease of finding and exploring 
new problems in numerical analysis is one of its particular attractions. 

As an illustration of this point that we hope will be of pedagogical interest to 
those teaching and thinking about numerical analysis, we pose and partially solve a 
question that was raised by the attractive cover picture of the book by Kincaid and 
Cheney (This book has been used at Duke University for the last few years in 
its graduate introductory course on numerical analysis; the first author was a recent 
student in this course when it was taught by the second author.) This cover picture 



is reproduced in Fig. 1.1 and reveals a remarkable result that is now nearly 120 years 
old jl], for the simple polynomial equation z 3 — 1 = in the complex plane, the 
sets of points zq that converge to a root under the map given by Newton's method, 
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have an unexpected intricate geometric complexity. In the language of nonlinear 
dynamics [ jl3| , the set of points that converge towards a fixed point under succes- 
sive iterations of some mapping is called the basin of attraction of that fixed point. 



Fig. 1.1 then says graphically that the basins of attraction for Newton's method are 
bounded by intricately braided objects. Naively, one would have expected the basins 
of attraction to be given by the sets of points closest to a given root and so bounded 
by straight lines. 



•This research was partially supported by grants from NSF NSF-DMS-93-07893 and DOE-DE- 
FG05-94ER25214 

tDepartment of Mechanical Engineering, Duke University, Box 90302, Durham, NC, 27708-0302 
bieOmel . egr .duke . edu 

tDepartment of Computer Science, Duke University, Box 90129, Durham, NC, 27708-0129 
hsgOcs . duke . edu 




Fig . 1.1. The regions of three different colors indicate the approximate basins of attraction for 
the three roots of the complex polynomial z [i — 1 = in the square region defined by — 1 < Real(z) < 1 
and — 1 < Imag(z) < 1. The positions of the roots of unity are indicated by the small black squares. 
This figure was computed numerically by using a mesh of 1000 X 1000 pixels in side this square. 
The center point of each pixel was used to start a Newton iteration under Eq. t \l.fy . The pixel 
was assigned a particular basin color if convergence to a root was identified within 1000 iterations. 
The criterion for convergence was that both the magnitudes of the correction 2;+i — Zi and of 
the residual z| — 1 were smaller than the value 10 — 12 . The boundaries separating the basins are 
intricately braided, with a fractal dimension of D m 1.42. 



This braided structure is typical of geometric objects known as fractals |2j, jL3|. 
The concept of a fractal is hard to define rigorously Q| but can be casually defined as 
a scale-invariant set that has a non-integer dimension D known as its fractal dimen- 



sion. For Fig. 1.1, scale invariance corresponds to the fact that if any small region 
of the braided structure is greatly magnified, further braided structure is found of a 
similar complexity and there is no end to this intricate detail upon further and further 
magnification. The fractal dimension is perhaps best understood as a scaling expo- 
nent that generalizes the usual concept of dimension to non-integer values |13|] . E.g., 
one kind of fractal dimension known as the capacity C indicates the scaling N e oc e~ c 
as e — > of the minimum number of balls N £ of radius e needed to cover the set as 
the ball radius e tends to zeroQ. Numerical calculations of the capacity of a set are 
often slowly convergent when the points of the set are visited with non-uniform prob- 
ability (measure) by the dynamics ML. For this reason, in the following we have used 
the average pointwise dimension D [|3| to quantify the structure of the fractal basin 
boundaries. This dimension takes into account the probability of visiting different 
parts of a set and so is more rapidly convergent in practice.^. 



1 More precisely, we would define the capacity as the limit C = lim e ^o log-/V(e)/ log(l/e) if this 
limit exists. In practice, one computes numerically a sequence of values N Ei for different radii Si 
and then estimates the dimension from the slope of a least-squares linear fit to a plot of log(7V e ) 
versus log(l/e). 

2 The point-wise dimension D is defined by the limit D = lim e ^o log-P(£)/ l°g(l/e) if this limit 
exists, where the quantity P(e) is the density or probability of points in a ball of radius e, i.e., the 
ratio of the number of points in the ball to the total number of points in the set. 



Fractal Basins of Attraction 



3 



Another interesting aspect of the fractal structure in Fig. [O] is that the bound- 
aries satisfy the counterintuitive Wada property []Io|| : every point on the boundary of 
one basin of attraction is also on the boundary of the other basins of attraction. 

The Kincaid and Cheney book touches only briefly on these remarkable facts and 
does not mention that fractal basins are in fact typical for Newton methods, rather 
than rare or pathological. Some of the important (and largely numerical) discoveries 
in the field of nonlinear dynamics in the last twenty years is that fractal structures 
pop up over and over again for nonlinear equations and for typical choices of equation 
parameters and of initial conditions: basin boundaries can be fractal, the attracting 
limit set for maps, flows, and partial differential equations can be fractal, and the 
set of parameter values that lead to a particular class of solutions (e.g., fixed points) 
can be fractal [fL3f . This fractal structure has important practical implications for 
many numerical problems and algorithms and so numerical analysis students should 
be aware of this. As one example, fractal basins of attraction for Newton's method 
imply that it will generally be difficult or impractical to characterize the ever-so- 
important initial guesses that will lead to convergence. In some extreme cases, it may 
be impossible to characterize the dependence of the algorithm or even of a laboratory 
experiment on parameters pT|. 



Fig. 1.1 is a good example of a result that is easily computed by students in 
an introductory numerical analysis course and that leads to interesting, and largely 
unexplored, questions that might entice a student to get interested in numerical anal- 
ysis. Thinking about this figure carefully raises the following kinds of questions: Are 
fractal boundaries associated with Newton's method rare or common? What de- 
termines the braided geometric structure (e.g., its fractal dimension D)l Do other 
root-finding algorithms lead to such fractal structure and, if so, can one relate the 
details of the algorithm to the fractal geometry? How do floating point errors affect 
the fractal structure? What are the implications of these fractal basin boundaries for 
root-finding software? 

In the following, we address briefly one of the above questions, namely do different 
root-finding algorithms — particularly Newton's method with damping — also lead to 
fractal basins of attraction? A well known weakness of Newton's method is that it 
is locally, but rarely globally, convergent Q. To improve this situation, researchers 
have invented damping methods |||5|] that try to increase the likelihood that an initial 
condition will lead to convergence. The basic strategy of most damping methods is to 
add only a fraction of a Newton step when correcting the current estimate of a root. 
During successive iterations, the damping algorithm usually increases the fraction 
of the Newton step towards the value of one until full Newton steps are taken and 
quadratic convergence is achieved. 

A question which, to our knowledge, seems not to have been discussed previously 
is then this: for fixed points with fractal basins of attraction, do damping methods 
change the fractal structure? If so, do they diminish or increase the fractal dimen- 
sion D of the boundaries? An interesting related question is whether damping works 
by increasing the basin of attraction, in that the basin with damping strictly contains 
the basin without damping. An increase in the fractal dimension D would possibly 
be an undesirable effect since more initial conditions would lie near the boundaries, 
making it harder to find criteria for identifying good initial conditions, which is often 
the most challenging part of a Newton method for many scientific problems. The 
question of how damping may modify the fractal basins of attraction is a natural one 
to ask given material normally covered in an introductory numerical analysis course 
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and is also within the ability of students to investigate themselves. 

2. The Fractal Basins Associated with Armijo's Rule. To study numeri- 
cally and graphically how the basins of attraction may be modified by damping, we 
consider one of the simplest damping methods, the Armijo rule For the equa- 

tion z 3 — 1 = 0, the Armijo rule modifies the Newton algorithm ( 1 . 1| ) by introducing 
a fractional damping coefficient 1/2 J in front of the Newton step as follows: 



(2.1) 
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where j is a nonnegative integer. Before the approximate root Zi is corrected to give a 
new and hopefully better estimate Zj+i, the integer j is increased from zero (i.e., the 
fraction of the step added is decreased) until the following criterion is just satisfied: 



(2.2) 
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At this point , Zi+i is computed from Eq. (2.1) with the given value of j and the Newton 
method is iterated until convergence or some maximum limit of iterations is attained. 
(If the latter, then usually a new initial condition must be tried.) Under certain 
assumptions of smoothness, one can prove that the Armijo method is guaranteed to 
converge to some root S . 
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Fig. 2 .1. Plot o f the basins of attraction of the roots of the equation z A — 1 = under Armijo's 

The fractal structure 
The basin boundaries 
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rule, Eqs. jj j and (2.1 ), calculated in the same way as described in Fig 
is preserved and slightly modified by damping when compared with Fig. 
have a fractal dimension D PS 1.29 with a standard deviation of 0.02 as determined from a least- 
squares fit to a plot o/log(P e ) versus log(l/e) (not shown). 



A comparison of the basins of attraction of the Newton- Armijo algorithm (Fig. |2.l| ) 
with those of the Newton algorithm (Fig. [Tl]) yields several interesting insights. First, 
the basin boundaries are not smoothed out under damping but retain a fractal struc- 
ture. Second, although the basins themselves are substantially altered (there are larger 
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contiguous regions of a given color in Fig. 2.1 compared to Fig. 1.1), damping causes 
only a modest change in the basin boundaries^ as quantified by the fractal dimen- 



sion D. The dimension of the fractal boundaries in Fig. 2.1 is D rj 1.29 whereas the 



fractal dimension in Fig. 1.1 has the somewhat larger value 1.42. The approximately 
ten percent relative decrease in the dimension seems consistent with the expected 
smoothing effect arising from the Armijo factor 1/2- 7 ' at each Newton- like step. Nev- 
ertheless, the dimension of the basin boundary does not generally decrease under 
damping as we point out in the next section. 

3. A Less Symmetric Example. Our discussion for the polynomial z 3 — 1 = 
is somewhat atypical in that there are only three roots and the corresponding basins 
of attraction are related by symmetry rotations around z = 0. To explore a more 
representative example and to suggest the rich geometric changes that damping can 
induce, we have also studied numerically the approximate basins of attraction of 
Newton and Newton- Armijo algorithms for the following two-dimensional nonlinear 
system: 

x — sin(x) cosh(w) = 0, 
(3.1) K ' y ' 

y — cos(x) sinh(y) = 0, 

which the second author first came across in a collection of numerical analysis prob- 
lems |l^, page 559, problem 16. E 



Figures 3.1 and 3.2 reveal a dramatic change in the shape and number of basins 
of attraction under damping although, within our numerical accuracy, the fractal 
dimension of a particular basin (the magenta basin belonging to the root (x, y) = 
(20.2, 3.7)) is unchanged with a value D » 1.4. 

Using starting points centered on 1500 x 1500 pixels spanning the region < 
x < 20 and < y < 5, we found that Newton's method converged^ towards 120 
different solutions^] Some of these solutions are shown in Fig. 3.3 and are color coded 



to indicate the associated basin of attraction. In contrast, the same initial conditions 
with Armijo's rule yielded only 53 different solutions, a substantially smaller number. 
These solutions were a subset of the 120 solutions obtained using Newton's method. 
Figures ET^EO raise fascinating open questions in addition to those already dis- 



cussed above in the context of Fig. 1.1. As an example, what has happened to the 
many basins of attraction that are now missing in the damped algorithm (53 solu- 
tions versus 120 for initial conditions in a certain given region)? Does the Wada 
property still apply to the various basins and, if so, which basins are linked by a 
given boundary? If not all basins are linked by the Wada property, do the different 
basin boundaries have different fractal dimensions and, if so, what determines these 
different dimensions? 



Finally, one can again ask how the geometry in Fig. 3.2 might depend on the 
damping algorithm. A canonical and interesting choice to study would be continuous 
damping in the sense of using an infinitesimal damping factor ds along a Newton 
direction. (We are indebted to our colleague, Dr. Donald Rose, for this suggestion.) 



3 A boundary was obtained by first separating the basin of attraction from the other basins, then 
by eliminating all the interior points, and finally, by computing the average pointwise dimension of 
the set of points that remained. 

4 The criterion for convergence was that the magnitude of the residual was smaller than 10 — 12 
and the magnitude of the correction term was smaller than 10 — 4 . 

5 Two solutions were considered different if the Euler distance between them in the x-y plane was 
greater than 10 -3 . 
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Fig. 3.1. Plot of the basins of attraction of the solutions of the system of equations (3.1) 
under Newton's method. Seven of the colors (green, purple, red, magenta, orange, yellow, and 
black) represent seven different basins of attraction (the most spatially extended regions). The blue 
regions represents all other basins while the white regions represent starting points for which the 
iterative process did not converge within 1000 iterations. The figure was calculated using a mesh 
of 1500 X 1500 pixels spanning the spatial region < x < 20 and < y < 5. The three solutions 
that lie inside this region are shown by the small black squares and have approximate coordinates 
(x,y) = (0,0), (7.5,2.8), and (13.9,3.4) respectively. The fractal dimension of the boundary of the 
magenta- colored basin (belonging to the approximate solution (x,y) = (20.2, 3. 7)) is D ss 1.38 ±0.04 
as determined from a least-squares fit to a plot o/log(P e ) versus log(l/e). 



Instead of iterating an Armijo-Newton map similar to Eqs. ( |2.l| ) and (2.2), one would 
then integrate a system of N ordinary differential equations given by: 

(3.2) « -j-i.ppq, 

as 

where the TV-dimensional system of equations is given by F(X) = 0, where J(X) is 
the N x TV Jacobian matrix J(X) = <9F/<9X, and where s is some arbitrary contin uou s 
parameter. Starting from some initial guess Xo, numerical integration of Eq. ( |3.2| ) 
will generate a path X(s) that must end up at a solution provided the Jacobian 
matrix remains nonsingular along this path. Challenging questions are raised in this 
continuous case: how long one must integrate to get close to a root? What integration 
accuracy is needed? What numerical algorithms should be used? What are the basins 
of attraction and are they fractal? 

4. Conclusions. Numerical analysis is an exciting and attractive subject specif- 
ically because it is fairly easy to to investigate novel problems even at the level of an 
introductory course. As an illustration of this, we have investigated numerically and 
graphically the problem of how a simple damping algorithm for Newton's method, 
Armijo's rule, modifies the fractal basins of attraction found for the complex poly- 




Fig. 3.2. The approximate basins of attraction of the solutions of Eq. ( [?. \ ) under Armijo's 
rule, using the same coloring scheme, resolution, and convergence criteria as in Fig. 3.1. The 
fractal dimension of the boundary of the magenta-colored basin colored is unchanged within numerical 
accuracy. 



nomial z 3 — 1 = . W e have also investigated a less symmetric and perhaps more 
general case, Eq. ( 3T ) , in which the basins of attraction are fundamentally changed 
upon damping. Our goal here was not to present a detailed or complete analysis but 
instead to suggest how students learning numerical analysis can easily find and explore 
new problems. The field of nonlinear dynamics [fl3f raises many related questions of 
geometric structure and their relation to numerical algorithms and is a valuable source 
of other ideas of pedagogical value for courses in numerical analysis. 
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